Near-infrared-laser-navigated dancing bubble within water via a thermally conductive interface

Precise manipulation of droplets or bubbles hosts a broad range of applications for microfluidic devices, drug delivery, and soft robotics. Generally the existing approaches via passively designing structured surfaces or actively applying external stimuli, inherently confine their motions within the planar or curved geometry at a slow speed. Consequently the realization of 3D manipulation, such as of the underwater bubbles, remains challenging. Here, during the near-infrared-laser impacting on water, by simply introducing a thermally conductive interface, we unexpectedly observe a spontaneously bouncing bubble with hundreds-of-micrometer diameter at tens-of-Hertz frequency. The unique formation of temperature inversion layer in our system generates the depth-dependent thermal Marangoni force responsible for the bouncing behavior. Both the scaling analysis and numerical simulation agree with observations quantitatively. Furthermore, by controlling the navigation speed of the laser beam, the bubble not only shows excellent steerability with velocity up to 40 mm/s, but also exhibits distinctive behaviors from bouncing to dancing within water. We demonstrate the potential applications by steering the bubble within water to specifically interact with tiny objects, shedding light on the fabrication of bubble-based compositions in materials science and contamination removal in water treatment.

laser beam and identify a range of bubble size (or of laser illu-mination duration) for which bouncing is observed. They provide measurements and simulations of the temperature and thermal convection flow fields induced by the laser, including in the thin temperature inversion layer below the heat-absorbent cover plate. They propose a bi-directional Marangoni-driven mechanism for the spontaneous bouncing of the bubble: downward Marangoni force due to the adverse temperature gradient below the plate and upward Marangoni force (adding up with buoyancy) due to the positive gradient in the bulk. They derive the scaling dependence of the lower and upper bounds of the bouncing bubble size range on the laser power, as well as the scaling dependence of the bouncing frequency on the bubble size, which they both compare to experimen-tal measurements.
Second, the authors present experiments with a translating laser beam, which displaces the bubble horizontally while preserving the vertical bouncing at low displacement speeds but not at large ones. They derive an expression for the maximal translating speed at which the bubble can follow the laser, as well as a scaling expression for the transition between bouncing and non-bouncing under translation. Finally the authors show examples where the displaced bubble can trap other bubbles, droplets or solid particles that float in the liquid, by coalescing, penetrating or wetting them.
The paper presents rather exhaustive observations and analyses of a generally interesting bouncing phenomenon as well as exemple of a potential use of it. However, although the different observations are gradually introduced and well related to each other, many crucial points in the presentation, in some experimental information and especially in the analyses and comparisons with the experimental data, which are listed below, are lacking or deserving signif-icant clarifications. Besides these crucial points, some of the discussions are difficult to follow, in particular technical discussions, because of syntactical issues in English (if the authors can clarify the crucial points, they should make sure to eliminate these syntactical issues). Last and importantly, the authors should clarify what is new and what is not among the many results they present.
1. Several arguments are not clear to me regarding the temperature evolution when the laser is not moving (S3).
I do not understand the scaling about the long-time temperature evolution (S3), ∆T ∝ t^(1/2). For a slender cylindrical laser beam with a steady power, the radial diffusion of heat should be 2 dimensional. In the absence of advection the temperature should reach a steady value, ∼ (r_l)^2*P/λ. Could the author comment and clarify which of the subsequent analyses this would modify and how? Also, could the authors clarify why the laser-induced convection can be neglected in this analyses?
2. The definition for the time origin (t = 0) is not explicitly given and it actually changes between the paper and the different supplementary materials without any explicit mention (between laser starting and bubble nucleation). This seems to be confusing not only for the reader but also for the authors themselves: The scaling for the laser induced temperature field (S3) cannot be used with a Plesset-Zwick law for the bubble growth without accounting for the time shift between the bubble nucleation and the laser starting. Could the author clarify?
3. Besides the two scalings mentioned above, the authors provide many scaling analyses but many of them are difficult or impossible to follow, because some notations are not introduced properly or not fully consistent, some expressions are not explicitly given (including in supplemental materials), and many scaling analyses only present the proportionality dependence on one parameter (e.g. laser power) while the other parameters are hidden.
a. The authors should define clearly each notation and give the actual experimental value for each of the parameters used in the analyses (for instance the value of dγ/dT, and all the other ones, should be given).
b. The full dimensional expression for each scaling law should be given before a specific scaling dependence on 1 parameter is highlighted.
c. For scaling dependences on 1 single parameter, the value of the dimensional pre-factor obtained from the theoretical analyses should be systematically given and compared to the experimentally observed value. Could the authors clarify how the temperature signal is quantitatively extracted from the IR camera signal and to which part of the liquid does this field actually refer to (given background thermal rays, reflexion, absorption and emissivity of the liquid)?

The derivation of the scaling for the thickness of the inversion layer is not clear (S7).
What is the situation really modeled for the convection velocity: a constant flux condition? applying on the surface of a cylinder? is the axial symmetry considered?
What is the dimensional factor missing in Vb ∝ Gr^(1/2)?
Also the final expression for δ_th is found to be independent on the plate properties, which seems contradictory with the authors' observation that the thermal effusivity of the plate is important. Could the authors clarify this point? 6. On page 9, it is not clear what is the status of the 0.1 in H_b,cr ∼ 0.1*R. Is this a fitting parameter or does it have a physical justification? In both cases this should be stated in the text. 7. On page 11 the derivation of eq (6) is impossible to follow. The assumptions should be clearly stated, the expression for τ_c and τ_d should be given, as well as all the necessary steps with all notations properly defined.
Also it would help the reader to summarize with simple arguments why and how the bouncing criteria differs when the laser moves relative to the case when it is static.
8. It seems obvious to the authors, though it is not explicitly justified, that the existence of a temperature inversion is sufficient to explain that the drop should bounce. However, one could also imagine that the bubble would reach a stable equilibrium position, close to the depth of the temperature maximum, where the net Marangoni force and buoyancy balance each other. It seems that the simulations presented in fig 3G and 4D do not answer to this point.
Could the authors comment? 9. Importantly, the authors should clarify what is new and what is not among the results they present. For instance, Marangoni induced displacement of bubbles by laser illumi-nation has already been reported (see e.g. Zhao et al, Lab Chip, 2014, 14, 384). Are equation (5) and the subsequent derivation of the maximal bubble displacement velocity original results? If not this should be explicitly stated. The originality should be clarified at each step in the presentation or analyses the authors report.
10. Relative to the demonstration in figure 5, how much is bouncing important in coalescence with bubbles, drops or particles? 11. In the introduction, I do not understand why Lohse & Zhang 2020 (ref 11) is cited for the bouncing of drops and bubbles. This review paper does not concern this thematics.
12. In the introduction, the seminal contribution regarding drops bouncing on a vibrated bath is the work of Couder and co-workers. It should, therefore, be cited.
13. The dimension is missing in the vertical axis in figure 4C.

Reviewer #3 (Remarks to the Author):
In this paper, the authors propose the use of thermal Maragoni flows to manipulate single bubbles and droplets. The results presented are compelling, the science is interesting, the research well executed, and the paper is rather well-written. Furthermore, the supplementary information brings real added value. I have, however, concerns which I detail below that I would like the authors to address before I can recommend publication in Nature Communications.
Motivation, claims and novelty: First, I have an issue with the motivation brought forth by the authors. In their listing of the technologies allowing for the manipulation of single bubbles or droplet, the authors mostly omit to discuss acoustical or optical tweezers which appear to be more straightforward, reliable, and impose less requirements on the sample. I do see the interest to create complex constructs as proposed in the third part of the paper. The claim for industrial relevance (p13 and introduction) also seems difficult to support, since the technique applies to the manipulation of single bubbles/droplets. Finally, oscillating bubbles induced by (thermal) Maragoni forces has been the object of several recent papers which is a concern for Nature Communications. Nevertheless, I do see the uniqueness of this particular physical configuration and I believe the novelty claims should be best focused on how this effect is put to use in the third part of the paper by further detailing the unique assemblies this technique can offer. other comments/concerns: "the observed 1/2 scaling relation of temperature evolution". I am puzzled by this statement.
First, based on the supplementary information (Fig. S3B), the scaling in ½ is not very convincing. In addition, ½ for heat diffusion (disregarding flow during transient build-up) corresponds to a plane, 1D cartesian diffusion. In Fig. S3B, we can see clearly the 1 power law for short time, corresponding to a (mostly) spherical diffusion of the heat deposited in the bulk by the gaussian Beam. In this context, isn't the decrease of slope observed after 1 second simply the slow transition toward the steady state temperature profile (slope = 0) existing for spherical heat diffusion? (Note that 1 s would be consistent with the typical diffusion timescale given the size of the laser beam). In this case, the effective power law would depend on the time at which nucleation occurs. If this is not the case, why does the power law in Fig S3B differs so significantly from ½ and why would the authors expect ½?. Finally, this scaling is bound to lose validity once the bubble is formed since it acts as a heat sink and changes the local absorption and light transmission. It is thus not trivial to me, based on the proposed arguments, why one would still expect a linear bubble growth over time. P4: "The existence of temperature inversion layer is attributed to the viscous and thermal boundary layers (BLs), as swept by the thermally-driven buoyancy flow during laser heating". What is dominant: the flow or the large conductivity of the sapphire that thus evacuates much heat near its surface? Thermal imaging: the authors do not provide much information on the thermal imaging. How do the authors solve the issue of imaging quantitively through a significant layer of water? In particular (p15), how is the thermal imaging calibrated? When the bubble is partly in and partly outside the boundary layer, can one simply consider 2 opposite forces, where the Maragoni force induced by the boundary layer is stronger (as suggested in the article) or does the bubble behaves as a dipole source for the flow? What would this difference entail? P2: "By utilizing an acoustic or magnetic field, non-invasive manipulation of bubbles is achieved, but sophisticated systems or potentially harmful sample pretreatments for additional properties are required" This is not necessarily true for acoustics-or optics-based techniques. Are the oscillations in Fig 1B real or only an effect of imperfect detection? If they are real, what is their source? P5: "the temperature peak is confirmed to be around hundreds of micrometers below the interface". Experimentally, the inversion goes until 0.4mm, this is twice as much as in the simulation. Furthermore, the profiles differ. Where do these differences come from? Finally, the experimental curve is plot in semilog and the simulation curve has linear axes. This should be uniform. Is the temperature profile decay for large z comparable in both cases? Fig 2E: PMMA has a lower conductivity than water. If, as the abstract suggests, thermal conductivity is key, why does it still induce a temperature inversion (albeit a small one)? Equation 4: although the scaling obtained seems to work, I am confused by the approach: the instability and thus the oscillation frequency in both cases originate from the competition of 2 opposing forces. Without either one, there would be no frequency. In term of oscillator, it is therefore the dependence on z of the resulting force that gives rise to the existence of an eigen frequency. How correct is then this analysis? This scaling implicitly considers both forces to be equal and opposite. I think these considerations should be more clearly introduced and discussed. P11: "the obtained maximum translation speed is consistent with the experimental observation, which is about 40 mm/s" How does this maximum translation speed depend on the laser power? P13: "the pre-existing bubble b and c along the pathway at 0.8 s, resulting in the sequential capture and fast coalescence at 0.9 and 1.0 s. Then the combined large bubble a + b + c continues its motion as guided by the laser (Movie S7)." I do not see the practical interest of this: there are many ways to induce coalescence, the difficulty is to prevent it. P14: "the solutal Marangoni force is vanished, and a distinct underlying mechanism for the bouncing motion is attributed to the thermal Marangoni" This statement is confusing: from the title of the reference discussed, it appears to also consider thermal Maragoni as (partial) driving mechanism. P16: section "PIV measurement for thermal buoyancy flow induced by laser heating". Since the PIV measurements are only shown in SI, I do not think it should be in the methods of the article itself. Fig 2D: the temperature profile does not seem to match Fig 2A. Why is that? The x-axes in Fig 2 B and D should have the same limits. Furthermore, these plots appear never to reach the state where the buoyancy force overcomes the Maragoni force. Eq. 2: where do these scaling laws exactly come from. P8, last line. The scaling for the thermal boundary layer thickness is obtained from this is a dimensional analysis with 3 nondimensional numbers. However, it entails non-trivial interconnected effects. Since the authors have thermal imaging capability, it would be interesting to verify this law experimentally. Furthermore, why would this law remain true in the presence of a bubble whose presence changes heat deposition and induces its own mixing? P9: "by setting Hb;cr = 0:1R." What is the physical argument behind this operation? Eq.4: why keeping the pi in the dimensional analysis when the 4/3 is already removed? Eq. 4: Equal sign of approximately equal?
Text: "followability": to me, this appear not to be the best choice of word here since the authors discuss the capacity of the bubble to follow the laser rather than the capacity of the laser to be followed. P2: 'fascinating'. I would suggest letting the reader judge of this. P3: "Here based on the technological advancement of laser impacting on water, by simply designing a specific thermally conductive interface, we directly realize both the vertical bouncing and the horizontal translating of the millimeter-sized bubble within pure water, achieving the dancing bubble within water in 3D." Please rephrase. P3: "Additionally, we demonstrate the remarkable manipulation capability of the bubble to interact with the preexisted droplet or nanoparticles" Please rephrase. P4: "During the whole operation process around 50 seconds," please rephrase P4: "robust or resilient" Do the authors mean "robust and resilient"? P4: "at t = t0" the plots refer to t=0. P6: "Such the temperature inversion might be responsible for the subsequent bouncing bubble." Please rephrase. P8: "two important issues are required to be addressed:" please rephrase

Reviewer 1
In the present work, authors observe a spontaneous bubble bouncing in a frequency of tens-of-Hertz by exposing water to a high power laser. To explain the underlying physical mechanism, authors conducted a systematic investigation through high-speed imaging, numerical simulation, PIV flow velocity measurement, scaling analysis of thermal Marangoni force. The bouncing behavior as well as the followability is verified by the simulation. Moreover, some of the potential applications are also demonstrated. Overall, the work is well implemented. I would recommend the acceptance of the manuscript after considering the following comments.
Thanks for "the acceptance of the manuscript after addressing the comments"! 1．As revealed by authors, the formation of temperature inversion layer is directly related to the high thermal conductivity of the sapphire glass. I wonder if the application of a top glass window with a low thermal conductivity will change the dynamics of bubble bouncing. Such a comparison will be very helpful.
Thanks for this wonderful comment! We have carefully elaborated the effect of the thermal conductivity of the cover materials by the following revision:  Clarify its physical mechanism in the section of "Formation of temperature inversion layer" in the main text.
In experiments, to further verify the effect of thermal conductivity of the top glass on the formation of temperature inversion layer, we tested various glasses with different thermal conductivity (see Table S1). The bubble bouncing behavior still holds for higher thermal conductivity cover such as quartz and sapphire cover, while bouncing disappears for lower thermal conductivity cover, such as PMMA and PDMS cover. As shown in Figure S3, no bouncing is observed for the PMMA cover, and bubble just hangs at the solid-liquid interface. In simulation, the corresponding simulated temperature distribution shown in Fig.  S6A also indicates that the thickness of temperature inversion layer becomes thinner for cover materials with a lower thermal effusivity, thus the temperature inversion layer is too thin to cause bubble bouncing for the case of PMMA and PDMS.
In theory, to qualitatively understand temperature inversion layer, we build the model to demonstrate the relevant physical parameters for the thickness of temperature inversion layer. The dependence of ratio between the thickness of temperature inversion layer ( ) and thickness of thermal boundary layer ( ℎ ) on the ratio between thermal effusivity of cover material ( ) and that of water ( ) is shown in Fig. S6C. The detail derivation is provided in the revised supplementary Section 7 "Model for temperature inversion layer". In the limit of high effusively or thermal conductivity (such as sapphire), the thickness of temperature inversion layer nearly approaches that of the thermal boundary layer ( ≈ ℎ ). Then, the thickness of temperature inversion layer extracted from simulation agree well with the scaling relation ∝ −1/4 for the analysis of the thickness of thermal boundary layer (Fig. S6D).
2．Water heats up by absorbing laser energy. The location of focal point along the vertical direction somehow will influence the generation of the temperature inversion layer and hence the dynamics of bubble bouncing. Can authors put some discussions on this?
Thanks for the comment. In our experiments, the laser beam output from a 980nm fiber coupled diode laser is first collimated by a plano-convex lens to a beam with diameter of about 2.2 mm. Then, the collimated laser beam is slightly focused by a focusing lens (L1) to a spot with diameter of about 1 mm.
The bouncing bubble is observed by either using the 2.2 mm collimated laser beam or the 1 mm slightly focused laser beam, although threshold laser intensity for bubble generation is different for the two cases.
Considering the penetration depth for 980-nm laser in water is about 2 cm, we treat the laser heating effect in our manuscript as a volumetric heating source, and assume that the location of focal point has a negligible effect on the formation of the temperature inversion layer (about 0.3 mm in our experiments).
In contrast, for a highly focused laser, the laser heating effect is treated as a point heat source, and the temperature peak is located at the focused point. The range of temperature distribution around bubble is about tens of micrometers, which is one order smaller than the thickness of temperature inversion layer in our experiments (about 0.3 mm), thus the resulting thermal Marangoni flow will trap the bubble to the focused spot [H. Takeuchi, Heat Transfer Engineering, 33, 234 (2012)], instead of bouncing.
In the revision, the sketch of laser in Fig. S1 is modified accordingly.   Thanks for the nice point! The symbol P indeed means the laser power with a unit of W in our manuscript, while laser intensity I is P/(πr 2 ) with a unit of W/m 2 . In the revised supplementary (Section 5), we correct the symbol P as I, which represents the laser intensity.
4．In the section "Mechanism of the bubble bouncing", authors used COMSOL to obtain numerical simulation results. They reconstructed the flow field. Figure 3G clearly shows that outside the thermal BLs, they obtained the Marangoni flow with the opposite direction. However, authors also claimed that outside the thermal BL, there is not temperature gradient. Can authors explain this inconsistency?
Thanks for the nice point! From the analysis in Fig. 3B in the main text and Section S8 in revised supplementary materials, when bubble is outside the temperature inversion layer, both the upward Marangoni force ( + ) and buoyancy force ( ) work as restoring forces. In order to reveal the flow patterns related with these driving forces, we perform simulations for three cases corresponding to (A) only + , (B) only , and (C) the combined two forces.
For case (A) only + , such as for a small bubble ( < ), since the upward Marangoni force + is much larger than buoyancy force , or the term of gravity in governing equation is ignored in analysis, and + dominates the rebouncing process. Correspondingly in simulation, the term of gravity in governing equation is ignored, and the positive temperature gradient outside the temperature inversion layer is assumed. The simulated flow pattern is obtained as shown in the following Figure  For case (B) only , such as for a large bubble ( > ), since buoyancy force is much larger than the upward Marangoni force + , + is ignored in analysis and dominates the rebouncing process. In simulation, the term of gravity in governing equation is included, while the temperature gradient outside the temperature inversion layer is neglected. The simulated flow pattern is obtained as shown in the following For case (C) the combined two forces, such as for a medium bubble ( ≈ ), both buoyancy force and the upward Marangoni force + should be considered. In simulation, we not only consider the term of gravity in governing equation, but also introduce a positive temperature gradient outside the temperature inversion layer. The simulated flow pattern is shown in the following figure (C), which is Fig. 3G in the revised main text.
The details of simulation settings and description for bubble bouncing are provided in the Section S11 in the revised supplementary materials.  Thanks for the nice comment! When some part of the bubble in the temperature inversion layer and the other part in the normal temperature gradient, there exists a competition of buoyancy force , upward Marangoni force + and downward Marangoni force − . Hence, in the revised Fig. 3D in the revised manuscript, we carefully compare the magnitude of these three forces: buoyancy force, upward thermal Marangoni force and downward Marangoni force. outside temperature inversion layer is about 1 K/mm. Thus, when the bubble touches surface, the downward Marangoni force is always at least one order larger than the upward Marangoni force. Therefore, to simplify the analysis, the upward Marangoni force can be reasonably neglected.
In the revised manuscript, we have presented the magnitude of the upward Marangoni force ( + ) in Fig. 3D to compare the individual component clearly. Thanks for the comment! Indeed, we also notice the lateral vibration at the early stage when the bubble radius is small.
A possible explanation for the lateral vibration is due to the Gaussian distribution of laser intensity, where the temperature peak along horizontal direction is located at center of the laser spot. The mass center of bubble tends to match the location of temperature peak, which is also the mechanism of bubble moving with the translating laser spot. However, for a stationary laser spot, bubble may be disturbed by the buoyancy flow of water, or the fluctuation of bubble radius, or the thermal-insulate effect of the bubble, to deviate from the center of the laser spot. Then, it appears to be lateral vibration.
In future, we will further study the lateral vibration with higher spatio-temporal resolution to characterize the exact motion of bubble at the early stage, and explore its underlying physical mechanism. In this work, we mainly focus on the periodic bouncing motion of bubble along the vertical direction. Figure  5, I did not see much relevance of the particle collection one to the manuscript. I recommend authors to remove this application.

7．Among the four applications of bubbles in composition fabrication shown in
Thanks for the comment! In literatures, bubble motion generally is passively manipulated along the translation direction of a laser beam, but in this work, we report the periodic bouncing motion of bubble along the penetration direction of laser beam. Besides the bouncing motion along the vertical direction, the floating bubble presented in this work can translate horizontally very fast with the laser spot (about 4 cm/s), which will enhance the efficiency and accuracy of particle collection.
Due to the extended moving dimension of bouncing bubble and the Marangoni convection around the bubble, the dancing bubble has a unique advantage in particle collection and transportation, hence in a fashion similar to "deep surface cleaning", particles can be collected or trapped by the bubble within a much larger volume (related with the vertical bouncing and horizontal translation). This novel capability of particle collection and transportation is demonstrated in the revised Supplementary Movie S10. Clearly, this remarkable particle enrichment and 3D manipulation of the bubble has implications for applications such as wastewater treatment or targeted drug delivery. Compared with other methods to trap particles, such as acoustically activated bubbles or previous optical manipulated bubbles, the underwater spherical dancing bubble in our system shows an obvious superiority in terms of enrichment concentration and flexible trajectory. Figure S2 should be depicted in Figure 1. Moreover, Figure S2B was not explained everywhere. The x axis scale in Figure S2B is not the same with that in Figure S2A.

8．Hb in
Thanks for the comment. In the revised manuscript, Hb is labeled as Ht , which is depicted in Fig. 1A. The description of Figure S2B is added in the revised supplementary (S2), also the x axis scale in Figure S2B is revised to be consistent with Figure S2A.
The following sentence is added in Supplementary Section S2 "Reproducibility of bouncing behavior" ---"The topmost position of the bubble Ht (as depicted in Fig. 1A, the vertical distance from the topmost of the bubble to the solid cover) is presented in Fig. S2B. The periodical displacement of the bubble implies the occurrence of the bouncing behavior, including the bouncing onset and stop moment, and the bouncing amplitude." 9．In the sixth line below figure 4 "the travel distance an be evaluated" "an"-> "can".
Thanks! This typo has been corrected, "the travel distance can be evaluated".

The paper presents experiments, scaling law arguments and numerical simulations of a phenomenon/process where a near-infrared laser beam is used to form, grow, bounce vertically, and displace horizontally a gas bubble in a liquid below a transparent top cover.
First, the authors present the bubble growth and bouncing phenomenology for a non-moving laser beam and identify a range of bubble size (or of laser illumination duration) for which bouncing is observed. They provide measurements and simulations of the temperature and thermal convection flow fields induced by the laser, including in the thin temperature inversion layer below the heat-absorbent cover plate. They propose a bi-directional Marangoni-driven mechanism for the spontaneous bouncing of the bubble: downward Marangoni force due to the adverse temperature gradient below the plate and upward Marangoni force (adding up with buoyancy) due to the positive gradient in the bulk. They derive the scaling dependence of the lower and upper bounds of the bouncing bubble size range on the laser power, as well as the scaling dependence of the bouncing frequency on the bubble size, which they both compare to experimental measurements.

Second, the authors present experiments with a translating laser beam, which displaces the bubble horizontally while preserving the vertical bouncing at low displacement speeds but not at large ones. They derive an expression for the maximal translating speed at which the bubble can follow the laser, as well as a scaling expression for the transition between bouncing and non-bouncing under translation.
Finally the authors show examples where the displaced bubble can trap other bubbles, droplets or solid particles that float in the liquid, by coalescing, penetrating or wetting them.

The paper presents rather exhaustive observations and analyses of a generally interesting bouncing phenomenon as well as example of a potential use of it.
Thanks for nice summary! We appreciate that "the paper presents rather exhaustive observations and analyses of a generally interesting bouncing phenomenon"! However, although the different observations are gradually introduced and well related to each other, many crucial points in the presentation, in some experimental information and especially in the analyses and comparisons with the experimental data, which are listed below, are lacking or deserving significant clarifications.
Besides these crucial points, some of the discussions are difficult to follow, in particular technical discussions, because of syntactical issues in English (if the authors can clarify the crucial points, they should make sure to eliminate these syntactical issues).

Last and importantly, the authors should clarify what is new and what is not among the many results they present.
Thanks for these constructive comments! These points will be addressed accordingly as below. Thanks for the comment. For a slender cylindrical laser beam with a steady power irradiating water, the temperature evolution has two limits: ∆ ∝ for short time period by ignoring thermal diffusion ( Fig. S4B and Fig. S4C), while ∆ ∝ 2 / for long time period by considering the balance of laser heating flux and thermal conductive flux. ) and the corresponding location (P=15 W with bubble formation) from the thermal images. (B, C) During the initial pre-heating stage, the evolution of elevated temperature (∆ ) for P=15 W with and P=10 W without bubble formation.
Actually, in previous supplementary materials, the long-time temperature evolution is used to explain the almost linear growth of bubble radius, but it does not affect other analysis for bubble motion. In the revision, in order to better focus on the bouncing behavior in this work, the analysis of 1/2 scaling of long-time temperature evolution and the linear growth of bubble are removed from the main text and supplementary materials in the revision, and the explanation of bubble motion is not affected.
By comparing the trivial solution for thermal conduction ∆ = 0 − 2 / 2 − (Eq. S4) with the simulation considering convection, as shown in Fig. R1, the temperature distribution outside the temperature inversion layer is well described by the trivial solution, while within the temperature inversion layer, the cooling effect of solid cover is dominant for our experiments. Hence, for the temperature evolution here, the laser-induced convection in the analysis for short time period can be ignored for simplicity.

The definition for the time origin (t = 0) is not explicitly given and it actually changes between the paper and the different supplementary materials without any explicit mention (between laser starting and bubble nucleation)
. This seems to be confusing not only for the reader but also for the authors themselves: The scaling for the laser induced temperature field (S3) cannot be used with a Plesset-Zwick law for the bubble growth without accounting for the time shift between the bubble nucleation and the laser starting. Could the author clarify?
Thanks for the nice comments! To avoid this unnecessary confusion, we present a consistent and clear definition of time origin in the revised supplementary materials. We define the time origin (t = 0) for the moment right before the bubble is visible, and the time for turning on laser during laser heating ( ), and we refine the related figures and descriptions in the revised main text and supplementary materials.
The different definition of time origin (the time shift between bubble nucleation and laser starting) can affect the scaling relation of temperature evolution. Besides the time shift effect, once the bubble forms, the thermal-insulation effect of bubble and the mixing effect of bubble motion also affect the temperature evolution as shown in Fig.  S4A after the time t = 0 s. Thus, we decide to remove the analysis of 1/2 scaling of long-time temperature evolution and the linear growth of bubble in the main text and supplementary materials. And the removed parts do not affect the explanation of bubble motion, which is the core part of this work. We will further study the temperature evolution and bubble growth characteristic in future. Thanks for the nice comment. In the revision, we have presented a clear definition for each notation, described each expression or parameters in detail, and added the value of parameters and their explicit and detailed expressions.
a. Definition and value for notations used in analysis is provided in Table S1 and  Table S2 the Supplementary.
b. The full dimensional derivation for each scaling laws are provided in the revised Supplementary Materials, including Section 6-9. Specifically, Section 6 is for "Model for velocity and thermal boundary layer"; Section 7 for "Model for temperature inversion layer"; Section 8 for "Model for bubble bouncing"; and Section 9 for "Model for the dancing bubble".
c. Comparison between theoretical estimation and experimental value for dimensional prefactors and physical parameters. Thanks for the nice comment! More details of the temperature measurements are provided here. In order to gather accurate temperature measurements, the temperature detection of the superheated liquid from the side view is calibrated, since the side wall of the quartz cuvette affects the infrared intensity signal to be detected by the camera from the target liquid surface.
Firstly, the radiometric thermal camera is so close to the target surface (less than 20 cm) that the atmospheric transmission factors can be neglected. Secondly, since the surface emissivity of both the quartz glass (ε=0.93) and liquid water (ε=0.96) is large (greater than 0.90), the impact of background temperature reflection can also be ignored. Thirdly, we calibrate the measured temperature from the IR camera by using a correction factor. The factor is obtained by simultaneously using the thermal camera and a thermocouple to measure the temperature evolution of a cooling process of hot glycerol (ε=0.96, the same emissivity as that of water) in the cuvette within the temperature range from 140 °C to 30 °C. And the location of the thermocouple is exactly the same with the laser focus spot, which is very near to the quartz wall (about several millimeter). Hence, the calibrated temperature refers to the part of the irradiated liquid in the plane of laser transmission.
A new section of "Temperature measurement" is added in Materials and Methods in the main text. "An infrared thermal camera (FLIR A6750sc) was used to acquire the temperature evolution of liquid during laser irradiating. For the top view temperature detection, the sapphire glass with high transmissivity @185-5000 nm was adopted as the cover glass since infrared light can optically penetrate the sapphire glass. Since the side wall of the quartz cuvette affects the detection of infrared intensity signal emitted from the target liquid surface, the measured temperature values from the side view by the thermal camera was required to be calibrated. For the calibration, a thermocouple was exactly located at the same position of laser focus spot, closer to the quartz wall (about several millimeter). By using the thermal camera and the thermocouple to simultaneously measure the temperature evolution of a cooling process for a glycerol in the cuvette within the temperature range from 140 °C to 30 °C (the same emissivity with that of water, ε = 0.96), the elevated temperature of the irradiated liquid in the plane of laser transmission was quantitatively calibrated." Thanks for this nice comment! We will clarify these points, including the effect of thermal conductivity of cover materials.

I.
For the model for the velocity and thermal boundary layer, detail derivation is provided in the revised supplementary Section S6 "Model for velocity and thermal boundary layer".
To  Clarify its physical mechanism in the section of "Formation of temperature inversion layer" in the main text.
In experiments, to further verify the effect of thermal conductivity of the top glass on the formation of temperature inversion layer, we tested various glasses with different thermal conductivity (see Table S1). The bubble bouncing behavior still holds for higher thermal conductivity cover such as quartz and sapphire cover, while bouncing disappears for lower thermal conductivity cover, such as PMMA and PDMS cover. As shown in Figure S3, no bouncing is observed for the PMMA cover, and bubble just hangs at the solid-liquid interface. In simulation, the corresponding simulated temperature distribution shown in Fig.  S6A also indicates that the thickness of temperature inversion layer becomes thinner for cover materials with a lower thermal effusivity, thus the temperature inversion layer is too thin to cause bubble bouncing for the case of PMMA and PDMS.
In theory, to qualitatively understand temperature inversion layer, we build the model to demonstrate the relevant physical parameters for the thickness of temperature inversion layer. The dependence of ratio between the thickness of temperature inversion layer ( ) and thickness of thermal boundary layer ( ℎ ) on the ratio between thermal effusivity of cover material ( ) and that of water ( ) is shown in Fig. S6C. The detail derivation is provided in the revised supplementary Section 7 "Model for temperature inversion layer".
In the limit of high effusively or thermal conductivity (such as sapphire), the thickness of temperature inversion layer nearly approaches that of the thermal boundary layer ( ≈ ℎ ). Then, the thickness of temperature inversion layer extracted from simulation agree well with the scaling relation ∝ −1/4 for the analysis of the thickness of thermal boundary layer (Fig. S6D). Thanks for the comment. As shown in Fig. 1C, we note that the bouncing begins with a small amplitude, to identify the lower threshold of bubble bouncing, a criterion for the bubble vertical displacement is needed. Thus, to eliminate the effect of fluctuation of bubble radius and lateral oscillation of bubble, we set the low bouncing threshold of , ≈ 0.1 that only when the bubble vertical displacement is comparable with 0.1 of its radius, based on the argument that bouncing becomes observable by high-speed images feasibly. The detail derivation is provided in the revised Supplementary Section 8. 7. On page 11 the derivation of eq (6) is impossible to follow. The assumptions should be clearly stated, the expression for τ_c and τ_d should be given, as well as all the necessary steps with all notations properly defined. Also it would help the reader to summarize with simple arguments why and how the bouncing criteria differs when the laser moves relative to the case when it is static.
Thanks for this comment! We have provided more details for equation (2) and equation (6) in the revised manuscript.
For equation (2) on page 8 in section of "Mechanism of the bubble bouncing", the related expressions are provided as below.
"When the bubble is outside the temperature inversion layer (Fig. 3A), both the buoyancy force = 4 3 3 ∝ 3 , and the upward Marangoni force + = ∆ • = ℎ + • 2 2 ∝ 2 ( ℎ ± = ± is the gradient of surface tension) serve as the restoring force, bringing bubble to the solid wall as expected. Once the bubble is mostly immersed in the temperature inversion layer (Fig. 3C), the buoyancy force is the same, and the upward Marangoni force + = ℎ + • max(0,2 − ) • depends on the height outside the temperature inversion layer. However, the downward Marangoni force − = ℎ − • min (2 , ) • becomes significant, enabling the possibility to push bubble downward away from the solid wall. From this scaling analysis, the comparison of the magnitude for these forces is shown in Fig. 3B  For equation (6) on page 11 in section of "Bubble dancing with horizontal translation", the related expressions are provided as below.
"Besides the dwelling timescale = 2 / , is the characteristic timescale for the heat transfer (defined in Fig. S4). For the short time period ( < ), the elevated temperature of water (∆ ) increases with time (∆ ∝ , Supplementary Section 5), while for the long time period ( > ), the elevated temperature of water reaches a steady state (∆ ∝ ).
Then, the following scaling is obtained, 3/4 / 2 ∝ [min ( , )] −1~{ −1~c onst.,( > , < 2 / ) −1~, ( < , > 2 / ) which agrees excellently with experiments (Fig. 4C) Thanks for this insightful comment! By fixing the temperature inversion layer in theory, bubble will finally tend to reach a stable equilibrium position. However, in our experiments, the formation of temperature inversion layer is associated with laser heating effect and convection flow, and might be disturbed by the thermal insulation effect of bubble and mixing effect induced by bubble bouncing. So, the temperature inversion layer actually experiences a construction-destruction-reconstruction process during bubble motion.
Although the laser is continuous, the downward Marangoni force impacted on bubble is in a fashion like a short impulse [similar idea to Zeng, et al, Proc. Natl. Acad. Sci. U. S. A. 118 (2021)]. In the simulation as shown in Fig 3G and 4D, the temperature inversion layer is introduced by a step function with time and only lasts for a short pulse (0< t <5 ms). By considering the stable temperature inversion layer and the short impulse acted temperature inversion layer, the simulated bubble motions are compared in the following figure. If the temperature inversion layer exists continuously or is fixed, the bubble reaches a stable equilibrium position in the bulk water.
Simulation details are provided in Sect 11 in the revised supplementary materials. Figure: The simulated bubble motion for two cases. (A) The temperature inversion layer acts as a short impulse on bubble motion. (B) The temperature inversion layer is fixed during the whole bubble motion process.

Importantly, the authors should clarify what is new and what is not among the results they present. For instance, Marangoni induced displacement of bubbles by laser illumination has already been reported (see e.g. Zhao et al, Lab Chip, 2014, 14, 384). Are equation (5) and the subsequent derivation of the maximal bubble displacement velocity original results? If not this should be explicitly stated. The originality should be clarified at each step in the presentation or analyses the authors report.
Thanks for the comment!
The main originality or novelty of this work is briefly summarized as following: (1) We observe the vertical bouncing bubble in water during laser impacting on water.
(2) We propose the underlying physical mechanism based on the unique temperature inversion layer, which is related with the high thermal conductivity of the cover glass. (3) We show the bouncing bubble can be steered horizontally by the laser, resulting in the dancing bubble. (4) We perform the scaling law analysis and numerical simulation, and the theory has a remarkable agreement with experiments.
In terms of Marangoni induced displacement of bubbles, previous work on bubble manipulation by laser is focused on bubble following movement with the translating laser spot or bubble trapping by a focused laser spot. Here, by designing the temperature field along laser propagation direction, we have successfully realized the bubble bouncing behavior, and this vertical bouncing behavior in water extends a new dimension along the vertical direction for bubble manipulation by laser.
In this work, we find bubble can bounce in water periodically along the direction of laser propagation, by only building a non-monotonic temperature distribution via laser heating. Moreover, by utilizing the horizontal motion of laser spot, we clearly demonstrate the 3D sophisticated manipulation of bubble or "dancing bubble" in water for the first time to the best of our knowledge.
As for the theoretical analysis, the full expression derivation and details are listed in the revised Supplementary Section 5-9, and the results referenced from previous work have been marked citations.

Relative to the demonstration in figure 5, how much is bouncing important in coalescence with bubbles, drops or particles?
Thanks for the comment! This work presents the bouncing bubble in water arising from the temperature inversion layer. This dancing bubble has two remarkable features, one is bouncing, and the other is steerability.
In the revised Fig 5, Fig5A and Fig 5B demonstrate these two features simultaneously, allowing 3D motion and leaping over a wall. Fig 5C-E mainly show its steerability to interact with the preexisting objects including walls, bubbles, droplets or nanoparticles, although bouncing might not be necessarily related.
Nevertheless, due to the extended moving dimension of bouncing bubble and the Marangoni convection around the bubble, the dancing bubble has a unique advantage in particle collection and transportation (Fig 5E) taken as an example, hence in a fashion similar to "deep surface cleaning", particles can be collected or trapped by the bubble within a much larger volume (related with the vertical bouncing and horizontal translation). This novel capability of particle collection and transportation is demonstrated in the revised Supplementary Movie S10. Clearly, this remarkable particle enrichment and 3D manipulation of the bubble has implications for applications such as wastewater treatment or targeted drug delivery. Compared with other methods to trap particles, such as acoustically activated bubbles or previous optical manipulated bubbles, the underwater spherical dancing bubble in our system shows an obvious superiority in terms of enrichment concentration and flexible trajectory.

In the introduction, I do not understand why Lohse & Zhang 2020 (ref 11) is cited
for the bouncing of drops and bubbles. This review paper does not concern this thematics.
Thanks for suggestions! This reference has been updated as below, "Particularly, arising from the effect of physicochemical hydrodynamics in a binary liquid (ethanol and water) (18)".

In the introduction, the seminal contribution regarding drops bouncing on a vibrated bath is the work of Couder and co-workers. It should, therefore, be cited.
Thanks for suggestion! The following two references have been cited accordingly---"such as droplets bouncing on a vibrating fluid bath (13-15)".

Review 3
In this paper, the authors propose the use of thermal Maragoni flows to manipulate single bubbles and droplets. The results presented are compelling, the science is interesting, the research well executed, and the paper is rather well-written. Furthermore, the supplementary information brings real added value.
Thanks for acknowledging "The results presented are compelling, the science is interesting, the research well executed, and the paper is rather well-written."

I have, however, concerns which I detail below that I would like the authors to address before I can recommend publication in Nature Communications.
Thanks for the insightful questions, which will be addressed as below.

Motivation, claims and novelty: First, I have an issue with the motivation brought forth by the authors. In their listing of the technologies allowing for the manipulation of single bubbles or droplet, the authors mostly omit to discuss acoustical or optical tweezers which appear to be more straightforward, reliable, and impose less requirements on the sample. I do see the interest to create complex constructs as proposed in the third part of the paper.
Thanks for this comment! This work mainly presents a physical phenomenon of bouncing bubble and dancing bubble during laser impacting on water underneath a glass cover, and proposes the underlying physical mechanism through scaling analysis and numerical simulation.
The comparisons with other technologies, acoustical and optical tweezers, will be discussed here. Firstly, optical tweezers have long been used as a powerful tool to trap high-index particles [A. Ashkin, Opt. Lett. 11, 288 (1986)], which rely on the extremely high gradient force produced near the beam waist of a tightly focused laser beam. However, a particle with a refractive index lower than its surrounding medium (a bubble being the limiting case) should be repulsed from the focus of a Gaussian beam. Although attempts have been made to overcome this limitation of conventional optical tweezers for trapping low-index particles by designing an optical vortex beam [K. T. Gahagan, Opt. Lett. 21, 827 (1996)] or using a two-dimensional interference pattern [M. P. MacDonald, Opt. Lett. 26, 863 (2001)], manipulation of bubbles in liquid by optical tweezers still remains challenging.
Secondly, acoustic methods are regarded as a promising alternative to trap and manipulate objects. Especially, the surface acoustic wave (SAW)-based microfluidic device has proven to be a potential powerful platform for manipulating bubbles, droplets and particles. In one-dimensional standing waves formed by two opposite travelling SAWs, the pressure gradient traps the bubbles at the potential wells, either the pressure nodes or pressure antinodes, according to the resonance frequency of bubbles with the respect to the acoustic frequency. Bubble position can be controlled by dynamically changing the position of the potential wells. Through designing a twodimensional standing surface acoustic waves with two pairs of transducers, the trajectories of bubbles can be configured in square or circular tracks [Long Meng, Appl. Phys. Lett. 100, 173701 (2012)]. However, technological limitations still exist for precise control of bubble location in three dimensions via these acoustic-based methods [Ding, Lab Chip, 2013, 13, 3626]. The manipulation in the out-of-plane (z) direction needs the innovation of device design and the improvement of the performance.
Additionally, a hybrid bubble-based manipulation technique called as "optoacoustic tweezers" for precise and on-demand handling of micro-objects by synergistically combining optothermal and acoustical strategies has been proposed [Xie, Lab Chip, 2013, 13, 1772. This technique utilizes the optothermal effect to generate size and location controllable bubbles and acoustic radiation forces to trap particles/cells.
Here, in our work, by only using the optothermal effect, the generation and 3D manipulation of bubbles, as well as particles trapping or concentrating are realized. Unlike optical or optoacoustic tweezers, our experimental setup does not require highnumerical aperture lenses or sophisticated systems. All setups are easy-to assemble and easy-to operate.
The comparison between these technologies has been added in the discussion part on page 14 in the revised main text, as below, "Conventionally, by utilizing the applied external stimuli (such as magnetic, electric, acoustic, and optical fields), the dynamics of droplets or bubbles can be accordingly controlled and their motions are subsequently directed (7-11). For example, magnetic and electric actuation require magneto-responsive and dielectric surfaces, but the need to manufacture patterning electrodes on substrate or embedding magnetic particles into a soft matrix increases the complexity of the overall operation (7, 8). Acoustic method offers the extraordinary capacity to trap and manipulate bubbles by forming the standing waves, but the precise control of bubble location in three dimensions still remains challenging (9, 10). Regarded as one of the promising methods for the remote and contactless manipulation, the optical approach such as the optical tweezers can trap high-refractive-index objects such as particles, but low-refractiveindex objects such as bubbles are anti-trapped (36)." The claim for industrial relevance (p13 and introduction) also seems difficult to support, since the technique applies to the manipulation of single bubbles/droplets.
Thanks for the comment! It is true that the current setup is difficult to support the industrial applications we referred to, but the proposed technique or method opens up a new possible way or sheds light on the relevant applications. For example, using bubbles for wastewater treatments or microorganism removal has been proven as a sustainable and environmentally friendly method in many industrial applications [T. Temesgen, Micro and nanobubble technologies as a new horizon for water-treatment techniques: A review, Adv. Colloid Interface Sci. 246, 40 (2017)].
Also, from a mechanics point of view, the shear stress when bubble impacting and sliding on a wall plays an important role in removing biofilms and contamination from various surfaces [Esmaili, Bubble impact on a tilted wall: Removing bacteria using bubbles, Phys. Rev. Fluid 4, 043603 (2019)]. Therefore, the laser manipulation of underwater bubble bouncing and moving is well suited for surface cleaning.
For material fabrication, the three-dimensional non-contact manipulation of underwater bubbles can induce bubbles to interact with other suspended or stationary bubbles/droplets for the construction of "gas-in-liquid" or "liquid-in-gas" composite capsules, such as the NIR-laser triggered core coalescence of double-emulsion drops.
Finally, since the Marangoni convection around thermal bubbles can drive the surrounding fluid to generate vortex patterns, which efficiently capture and enrich microparticles on the bubble surface for carrying and transferring drug molecules, cells, microorganisms. This particle enrichment and 3D manipulation of the bubble might be relevant for applications such as wastewater treatment or targeted drug delivery.
Actually, the technique we propose here is not only applicable for manipulation of single bubble/droplet, but also possibly feasible for multiple bubbles. The single laser beam can be easily split into several beams or programmed using a beam splitting element (diffractive optical elements) or spatial light modulators. Thus, generation and manipulation of multi-bubbles simultaneously can be achieved, likely holding promising for the industrial application.
Finally, oscillating bubbles induced by (thermal) Maragoni forces has been the object of several recent papers which is a concern for Nature Communications. Nevertheless, I do see the uniqueness of this particular physical configuration and I believe the novelty claims should be best focused on how this effect is put to use in the third part of the paper by further detailing the unique assemblies this technique can offer.

Thanks for the comment!
The main originality or novelty of this work is briefly summarized as following: (1) We observe the vertical bouncing bubble during laser impacting on water.
(2) We propose the underlying physical mechanism based on the unique temperature inversion layer, which is related with the high thermal conductivity of the cover glass. (3) We show the bouncing bubble can be steered horizontally by the laser, resulting in the dancing bubble. (4) We perform the scaling law analysis and numerical simulation, and the theory has a remarkable agreement with experiments.
In terms of the potential applications, previous work on bubble manipulation by laser is focused on bubble following the moving laser spot or bubble trapped by a focused laser spot. Here, by designing the temperature field along laser propagation direction, we have successfully realized predictable bubble bouncing behavior consequently extending a new dimension (along the propagation direction) for bubble manipulation by laser. This dancing bubble in our work has two remarkable features, one is bouncing, and the other is steerability (Fig .5). Fig5A and Fig 5B demonstrate these two features simultaneously, allowing 3D motion and leaping over a wall. Fig 5C-E mainly show its steerability to interact with the preexisting objects including walls, bubbles, droplets or nanoparticles, although bouncing might not be necessarily related.
Nevertheless, due to the extended moving dimension of bouncing bubble and the Marangoni convection around the bubble, the dancing bubble has a unique advantage in particle collection and transportation (Fig 5E) taken as an example, hence in a fashion similar to "deep surface cleaning", particles can be collected or trapped by the bubble within a much larger volume (related with the vertical bouncing and horizontal translation). This novel capability of particle collection and transportation is demonstrated in the revised Supplementary Movie S10. Clearly, this remarkable particle enrichment and 3D manipulation of the bubble has implications for applications such as wastewater treatment or targeted drug delivery. Compared with other methods to trap particles, such as acoustically activated bubbles or previous optical manipulated bubbles, the underwater spherical dancing bubble in our system shows an obvious superiority in terms of enrichment concentration and flexible trajectory.
other comments/concerns: 1. "the observed 1/2 scaling relation of temperature evolution". I am puzzled by this statement. First, based on the supplementary information (Fig. S3B), the scaling in ½ is not very convincing. In addition, ½ for heat diffusion (disregarding flow during transient build-up) corresponds to a plane, 1D cartesian diffusion. In Fig. S3B, we can see clearly the 1 power law for short time, corresponding to a (mostly) spherical diffusion of the heat deposited in the bulk by the gaussian Beam. In this context, isn't the decrease of slope observed after 1 second simply the slow transition toward the steady state temperature profile (slope = 0) existing for spherical heat diffusion? (Note that 1 s would be consistent with the typical diffusion timescale given the size of the laser beam). In this case, the effective power law would depend on the time at which nucleation occurs.
If this is not the case, why does the power law in Fig S3B differs so significantly from ½ and why would the authors expect ½?.
Finally, this scaling is bound to lose validity once the bubble is formed since it acts as a heat sink and changes the local absorption and light transmission. It is thus not trivial to me, based on the proposed arguments, why one would still expect a linear bubble growth over time.
Thanks for the comment. For a slender cylindrical laser beam with a steady power irradiating water, the temperature evolution has two limits: ∆ ∝ for short time period by ignoring thermal diffusion ( Fig. S4B and Fig. S4C), while ∆ ∝ 2 / for long time period by considering the balance of laser heating flux and thermal conductive flux. Actually, in previous supplementary materials, the long-time temperature evolution is used to explain the almost linear growth of bubble radius, but it does not affect other analysis for bubble motion. In the revision, in order to better focus on the bouncing behavior in this work, the analysis of 1/2 scaling of long-time temperature evolution and the linear growth of bubble are removed from the main text and supplementary materials in the revision, and the explanation of bubble motion is not affected.

P4:"The existence of temperature inversion layer is attributed to the viscous and thermal boundary layers (BLs), as swept by the thermally-driven buoyancy flow during laser heating". What is dominant: the flow or the large conductivity of the sapphire that thus evacuates much heat near its surface?
Thanks for this wonderful comment! We have clarified the underlying physical mechanism of temperature inversion layer in the revision.
In short, thermal boundary layer ( ℎ ) is caused by the thermally-driven buoyancy flow. The temperature inversion layer indicates the temperature distribution is nonmonotonic, and its thickness ( ) is determined by determined by the competition between laser heating rate and cover cooling rate. More quantitatively, / ℎ depends on the ratio of their thermal effusivity between cover materials and water, as shown below (the same as Figure S6C in Supplementary). Figure S6C: The thickness of temperature inversion layer dependent on the thermal effusivity of cover material Below we will clarify these points, particularly the effect of thermal conductivity of cover materials on .
I. For the model for the velocity and thermal boundary layer ( ℎ ), detail derivation is provided in the revised supplementary Section S6 "Model for velocity and thermal boundary layer".
To model the boundary layer induced by the buoyancy flow sweeping the cover surface, we first consider the buoyancy flow similar to the analysis for laminar natural convection on a heated vertical surface [Fundamentals of Heat and Mass Transfer (John Wiley & Sons, 2011)], which is an axisymmetric model. The expression of buoyancy flow velocity is obtained = √ ∆ = 1/2 . Then, we use the hydrodynamic solution by Blasius for laminar flow over isothermal plate [Fundamentals of Heat and Mass Transfer (John Wiley & Sons, 2011)] to obtain the expression of velocity and thermal boundary layer thickness.

II.
We have carefully elaborated the effect of the thermal conductivity of the cover materials by the following revisions:  Perform experiments for various cover materials, as shown in Supplementary Section 3 "The effect of cover materials on the bouncing behavior" with Fig. S3;  Build model to understand temperature inversion layer, as shown in Supplementary Section 7 "Model for temperature inversion layer" with Fig. S6;  Clarify its physical mechanism in the section of "Formation of temperature inversion layer" in the main text. In experiments, to further verify the effect of thermal conductivity of the top glass on the formation of temperature inversion layer, we tested various glasses with different thermal conductivity (see Table S1). The bubble bouncing behavior still holds for higher thermal conductivity cover such as quartz and sapphire cover, while bouncing disappears for lower thermal conductivity cover, such as PMMA and PDMS cover. As shown in Figure S3, no bouncing is observed for the PMMA cover, and bubble just hangs at the solid-liquid interface.
In simulation, the corresponding simulated temperature distribution shown in Fig.  S6A also indicates that the thickness of temperature inversion layer becomes thinner for cover materials with a lower thermal effusivity, thus the temperature inversion layer is too thin to cause bubble bouncing for the case of PMMA and PDMS.
In theory, to qualitatively understand temperature inversion layer, we build the model to demonstrate the relevant physical parameters for the thickness of temperature inversion layer. The dependence of ratio between the thickness of temperature inversion layer ( ) and thickness of thermal boundary layer ( ℎ ) on the ratio between thermal effusivity of cover material ( ) and that of water ( ) is shown in Fig. S6C. The detail derivation is provided in the revised supplementary Section 7 "Model for temperature inversion layer".
In the limit of high effusively or thermal conductivity (such as sapphire), the thickness of temperature inversion layer nearly approaches that of the thermal boundary layer ( ≈ ℎ ). Then, the thickness of temperature inversion layer extracted from simulation agree well with the scaling relation ∝ −1/4 for the analysis of the thickness of thermal boundary layer (Fig. S6D).

Thermal imaging: the authors do not provide much information on the thermal imaging. How do the authors solve the issue of imaging quantitatively through a significant layer of water? In particular (p15), how is the thermal imaging calibrated?
Thanks for the nice comment! More details of the temperature measurements are provided here. In order to gather accurate temperature measurements, the temperature detection of the superheated liquid from the side view is calibrated, since the side wall of the quartz cuvette affects the infrared intensity signal to be detected by the camera from the target liquid surface.
Firstly, the radiometric thermal camera is so close to the target surface (less than 20 cm) that the atmospheric transmission factors can be neglected. Secondly, since the surface emissivity of both the quartz glass (ε=0.93) and liquid water (ε=0.96) is large (greater than 0.90), the impact of background temperature reflection can also be ignored. Thirdly, we calibrate the measured temperature from the IR camera by using a correction factor. The factor is obtained by simultaneously using the thermal camera and a thermocouple to measure the temperature evolution of a cooling process of hot glycerol (ε=0.96, the same emissivity as that of water) in the cuvette within the temperature range from 140 °C to 30 °C. And the location of the thermocouple is exactly the same with the laser focus spot, which is very near to the quartz wall (about several millimeter). Hence, the calibrated temperature refers to the part of the irradiated liquid in the plane of laser transmission.
A new section of "Temperature measurement" is added in Materials and Methods in the main text. "An infrared thermal camera (FLIR A6750sc) was used to acquire the temperature evolution of liquid during laser irradiating. For the top view temperature detection, the sapphire glass with high transmissivity @185-5000 nm was adopted as the cover glass since infrared light can optically penetrate the sapphire glass. Since the side wall of the quartz cuvette affects the detection of infrared intensity signal emitted from the target liquid surface, the measured temperature values from the side view by the thermal camera was required to be calibrated. For the calibration, a thermocouple was exactly located at the same position of laser focus spot, closer to the quartz wall (about several millimeter). By using the thermal camera and the thermocouple to simultaneously measure the temperature evolution of a cooling process for a glycerol in the cuvette within the temperature range from 140 °C to 30 °C (the same emissivity with that of water, ε = 0.96), the elevated temperature of the irradiated liquid in the plane of laser transmission was quantitatively calibrated." 4. When the bubble is partly in and partly outside the boundary layer, can one simply consider 2 opposite forces, where the Maragoni force induced by the boundary layer is stronger (as suggested in the article) or does the bubble behaves as a dipole source for the flow? What would this difference entail?
Thanks for the nice comment! When some part of the bubble in the temperature inversion layer and the other part in the normal temperature gradient, there exists a competition of buoyancy force , upward Marangoni force + and downward Marangoni force − . Hence, in the revised Fig. 3D in the revised manuscript, we carefully compare the magnitude of these three forces: buoyancy force, upward thermal Marangoni force and downward Marangoni force.
For the case that bubble touches surface (as shown in Fig. 3C outside temperature inversion layer is about 1 K/mm. Thus, when the bubble touches surface, the downward Marangoni force is always at least one order larger than the upward Marangoni force. Therefore, to simplify the analysis, the upward Marangoni force can be reasonably neglected. In the revised manuscript, we have presented the magnitude of the upward Marangoni force ( + ) in Fig. 3D to compare the individual component clearly. and their magnitudes dependent on bubble radius 5. P2: "By utilizing an acoustic or magnetic field, non-invasive manipulation of bubbles is achieved, but sophisticated systems or potentially harmful sample pretreatments for additional properties are required" This is not necessarily true for acoustics-or optics-based techniques.
Thanks for the comment! We have refined the expression in Introduction on page 2, "Other possible approaches utilize the external stimuli to actively drive the motion of droplets or bubbles, including magnetic, electric, acoustic or light fields (7-11), and their dynamics is correspondingly restricted by the applied specific forces".
We have included the comparison of these technologies in Discussion on page 14, "Conventionally, by utilizing the applied external stimuli (such as magnetic, electric, acoustic, and optical fields), the dynamics of droplets or bubbles can be accordingly controlled and their motions are subsequently directed (7-11). For example, magnetic and electric actuation require magneto-responsive and dielectric surfaces, but the need to manufacture patterning electrodes on substrate or embedding magnetic particles into a soft matrix increases the complexity of the overall operation (7,8). Acoustic method offers the extraordinary capacity to trap and manipulate bubbles by forming the standing waves, but the precise control of bubble location in three dimensions still remains challenging (9, 10). Regarded as one of the promising methods for the remote and contactless manipulation, the optical approach such as the optical tweezers can trap high-refractive-index objects such as particles, but low-refractive-index objects such as bubbles are anti-trapped (36)." 6. Are the oscillations in Fig 1B real  Thanks for the comment! The oscillations of bubble radius in Fig. 1B can be divided into two kinds: one is the increased steps, which is caused by the limitation of spatial resolution, the other is strong oscillations during bubble bouncing, which may be caused by the bubble deformation and the inner pressure fluctuation varied with the surrounding liquid temperature. We identify the bubble bouncing based on the timedependent position of bubble center (Hc as shown in Fig. 1C). 7. P5: "the temperature peak is confirmed to be around hundreds of micrometers below the interface". Experimentally, the inversion goes until 0.4mm, this is twice as much as in the simulation. Furthermore, the profiles differ. Where do these differences come from? Finally, the experimental curve is plot in semilog and the simulation curve has linear axes. This should be uniform. Is the temperature profile decay for large z comparable in both cases?
Thanks for the comment! In the revised manuscript, we uniform the x axis ( Fig.  2B and E) in linear scale, and the simulated temperature data for sapphire with thickness of 1 mm (blue line in previous Fig. 2E) is replaced by data for sapphire with thickness of 250 μm (blue line in revised Fig. 2E) which is consistent with the glass thickness (250 μm) in experiment.
The thickness of temperature inversion layer ( ) is fluctuating during the bouncing process but has an average value of about 0.3 mm, as shown in Fig. S3A. In the revised Fig. 2B and E, the thickness of temperature inversion layer is about 0.3 mm experimentally and numerically, also their profiles are similar.
Since temperature profile is simulated to validate the existence of the temperature inversion layer due to the high thermal conductivity of the cover glass, temperature profile decays for large Z is not focused here.
8. Fig 2E:  Thanks for the comment. The bouncing bubble system is more like the parabolic motion rather than the oscillator. The downward Marangoni force induced by temperature inversion layer provides the initial velocity of bubble motion. Then, the motion of bubble destructs the temperature inversion layer. The buoyancy force and upward Marangoni force (induced by the Beer-Lambert law) act as the restoring forces, then the bubble returns back to the cover, and ends a cycle of bounce.
Thus, the oscillation frequency in our experiments is obtained not from eigen frequency, but from solving the bubble motion equation. This idea based on the bubble motion equation is similar to study bouncing bubble induced by the competition between driving solute Marangoni force and restoring thermal Marangoni force in a recent literature [Zeng, et al, Proc. Natl. Acad. Sci. U. S. A. 118 (2021)]. Furthermore, the detail derivation process is provided in the revised Supplementary Section 8.
Also we compare the estimated value in analysis and the observed value in experiments. For small bubble, the upward Marangoni force + is dominant, and the obtained expression of frequency = 1 −1/2 . The estimated value of the prefactor 1 ≈ 10, and the observed value of the prefactor 1 ≈ 16 in Fig. 3F. For large bubble, the buoyancy force is dominant. For bubble radius = 0.9 mm, the estimated frequency = 19.6 Hz, and the observed frequency is = 19.2 ± 2.17 Hz in Fig. 3F To prevent bubble coalescence, our method to obtain dancing bubble might be relevant by jumping over the pre-existing bubble, but requires more future work. 12. P14: "the solutal Marangoni force is vanished, and a distinct underlying mechanism for the bouncing motion is attributed to the thermal Marangoni" This statement is confusing: from the title of the reference discussed, it appears to also consider thermal Maragoni as (partial) driving mechanism.
Thanks for the comment! In this reference, they use solute Marangoni force as a driving force, and thermal Marangoni force as a restoring force, which is different from the situation in our manuscript that using bidirectional thermal Marangoni force and buoyancy force to achieve bouncing.
To clarify the unnecessary confusion, the sentence has been modified in the main text on page 15 as below, "We note that in the recent work in Ref (23), bouncing plasmonic bubble has been demonstrated in a binary liquid consisting of water and ethanol, and the competition between the solutal and thermal Marangoni forces is identified as the origin of the periodic bouncing. However, with only pure water as the host fluid here, both the bidirectional thermal Marangoni force and buoyancy force are implemented carefully to achieve bouncing.". Thanks for suggestion! This part has been removed in the methods and included in the Supplementary Section 6 in the revised manuscript.
14. Fig 2D: the temperature profile does not seem to match Fig 2A. Why is that? Thanks for the comment. As limited by the spatial resolution of thermal camera, the snapshots from thermal camera is used to show the overall distribution. The simulation snapshot is a zoom-in one to show more detail within temperature inversion layer. The differences in experimental and numerical temperature profile might be due to the effect of bubble formation and bounce. Fig 2 B and D should have the same limits. Furthermore, these plots appear never to reach the state where the buoyancy force overcomes the Maragoni force.

The x-axes in
Thanks for the comment! We have updated the corresponding figure accordingly ( Fig. 3B and D). When the bubble radius exceeds , the buoyancy force overcomes the downward Marangoni force − , and the bubble cannot bounce any more. Thanks for the nice comment! When some part of the bubble in the temperature inversion layer and the other part in the normal temperature gradient, there exists a competition of buoyancy force , upward Marangoni force + and downward Marangoni force − . Hence, in the revised Fig. 3D in the revised manuscript, we carefully compare the magnitude of these three forces: buoyancy force, upward thermal Marangoni force and downward Marangoni force.
For the case that bubble touches surface (as shown in Fig. 3C outside temperature inversion layer is about 1 K/mm. Thus, when the bubble touches surface, the downward Marangoni force is always at least one order larger than the upward Marangoni force. Therefore, to simplify the analysis, the upward Marangoni force can be reasonably neglected. In the revised manuscript, we have presented the magnitude of the upward Marangoni force ( + ) in Fig. 3D to compare the individual component clearly. We have included the laser power in the caption of Fig. 3 (P = 20 W) . Fig 3 B and D.

It would be helpful if Rup and Rlow are depicted in
Thanks! We have depicted Rup in Fig. 3D, but Rlow is not from the force analysis.
Eq. 2: where do these scaling laws exactly come from.
In the revised manuscript, the full expressions of + and − are given. The P8, last line. The scaling for the thermal boundary layer thickness is obtained from this is a dimensional analysis with 3 nondimensional numbers. However, it entails nontrivial interconnected effects. Since the authors have thermal imaging capability, it would be interesting to verify this law experimentally.
Furthermore, why would this law remain true in the presence of a bubble whose presence changes heat deposition and induces its own mixing?
Details of the scaling for the thermal boundary layer thickness can be seen in the revised supplementary section 6 ("Model for velocity and thermal boundary layer").
To model the boundary layer induced by the buoyancy flow sweeping the cover surface, we first consider the buoyancy flow similar to the analysis for laminar natural convection on a heated vertical surface [Fundamentals of Heat and Mass Transfer (John Wiley & Sons, 2011)], which is an axisymmetric model. The expression of buoyancy flow velocity is obtained = √ ∆ = 1/2 . Then, we use the hydrodynamic solution by Blasius for laminar flow over isothermal plate [Fundamentals of Heat and Mass Transfer (John Wiley & Sons, 2011)] to obtain the expression of velocity and thermal boundary layer thickness. Then, the thickness of velocity boundary layer = 5 /√ , and the thickness of thermal boundary layer ℎ = / 1/3 . The physical relation between the thickness of temperature inversion layer ( ) and that of thermal boundary layer ( ℎ ) has been clarified in the revised Supplementary Section 7.
In fact, the formation of temperature inversion layer in our experiments is associated with laser heating effect and stable convection fluid flow (about 0.01 m/s in experiments), while the flow induced by bubble motion (~ 0.1 m/s) is one order larger than the buoyancy flow. Then, the temperature inversion layer actually experiences a construction-destruction-reconstruction process during bubble motion. This complicated dynamical process hinders to verify the formula experimentally.
Thanks for the comment. As shown in Fig. 1C, we note that the bouncing begins with a small amplitude, to identify the lower threshold of bubble bouncing, a criterion for the bubble vertical displacement is needed. Thus, to eliminate the effect of fluctuation of bubble radius and lateral oscillation of bubble, we set the low bouncing threshold of , ≈ 0.1 that only when the bubble vertical displacement is comparable with 0.1 of its radius, based on the argument that bouncing becomes observable by high-speed images feasibly. The detail derivation is provided in the revised Supplementary Section 8.

Eq.4: why keeping the pi in the dimensional analysis when the 4/3 is already removed?
We have removed this coefficient accordingly.
Text: "followability": to me, this appear not to be the best choice of word here since the authors discuss the capacity of the bubble to follow the laser rather than the capacity of the laser to be followed.
The word "followability" is replaced by "steerability" in the revised manuscript. We have removed the word "fascinating".
P3: "Here based on the technological advancement of laser impacting on water, by simply designing a specific thermally conductive interface, we directly realize both the vertical bouncing and the horizontal translating of the millimeter-sized bubble within pure water, achieving the dancing bubble within water in 3D." Please rephrase.
The sentence has been rephrased as below "Here, by adopting a near-infrared laser irradiating into pure water through a transparent glass cover with a high thermal conductivity, we observe a vertical bouncing behavior of an optothermal bubble…Moreover, being subject to the navigation of the laser beam, the resulted dancing bubble can be translated along the horizontal direction." P3: "Additionally, we demonstrate the remarkable manipulation capability of the bubble to interact with the preexisted droplet or nanoparticles" Please rephrase.
The sentence has been rephrased as below "Further, through the precise control of the specific interaction with the preexisting objects, we successfully demonstrate the unique bouncing and steerability features of this dancing bubble".
P4: "During the whole operation process around 50 seconds," please rephrase The sentence has been rephrased as below "During the bubble growth within 50 s" P4: "robust or resilient" Do the authors mean "robust and resilient"?
This sentence has been removed in the revised manuscript.
We define t=0 as the moment right before the bubble is visible, and t= t0 = 18.48 s in Fig. 2A and 2B refers to the onset moment of one bouncing cycle. The sentence " Fig.  2A and B at t=t0 and t0 + 56 ms" has been removed in the revision accordingly.
P6: "Such the temperature inversion might be responsible for the subsequent bouncing bubble." Please rephrase.
This sentence has been rephrased as below, "Particularly, for the glass cover with a high thermal conductivity, within TIL, the temperature gradient up to 10 K/mm (Fig.  S4) might be responsible for the bouncing behavior." P8: "two important issues are required to be addressed:" please rephrase Thanks! This sentence has been removed in the revised main text.

P9: "converted into kinetic energy of bubble with the initial velocity" please rephrase
This sentence has been rephrased as below "From energy perspective, the work done by the downward Marangoni force ( − ∝ ℎ − 2 ) is mainly converted into the kinetic energy of bubble associated with an initial velocity ( ,0 ), ∝ I thank the authors for their detailed answers and the many clarifications they have brought to both the paper and the supplementary materials. They have clarified the major concerns I had. Their analysis and discussions are, now, presented in a clear, refutable and convincing way. I think the work they present is worth being published.
Thanks for acknowledging the improvement of the manuscript, and we appreciate that "the work they present is worth being published"! I think they should still adjust one explanation regarding point 1 of my previous report, which does not change the following of their analysis. The bouncing phenomenon is observed for the 'long-time' regime of the heating of the liquid by the laser. Therefore, equations S3 and S4 are not relevant. They should rather be replaced by the steady, radial heat equation (1/r)d(rκ dT /dr)/dr + A0P exp(−r2/rl2) = 0 and its solution, which is also linear in P.
Thanks for the nice point! The analysis of steady temperature for the "long-time" regime is added in Supplementary Discussion 1 "Model for peak temperature evolution" in the revised supplementary materials as suggested.
"For temperature evolution during the long time period, by neglecting convection and considering the thermal diffusion along the r direction, the steady-state axisymmetric heat conduction equation is, 1 ( ) + 0 − 2 / 2 = 0 The boundary conditions are as below, Here the first boundary condition results from the symmetry, and the second condition is a constant temperature for the far field condition, where rw is the width of water bulk (rw= 5 mm in experiments).
For sufficiently long time period, the diffusion length of heat transfer is much longer than the laser beam, and the laser intensity profile (the term − 2 / 2 ) is simplified to a point source of heat. Then, according to Ref. (2), the analytical solution is, This also suggests that the authors should simply remove all detailed discussions regarding the short-time (unsteady) regime when the liquid temperature is increasing. This is unnecessary and may bring confusion.
Thanks for the nice comment! Since the short-time (unsteady) regime is responsible for the analysis for bubble dancing when laser is translated horizontally, as shown in Supplementary Discussion 5 "Model for the dancing bubble". Hence, in the revised supplementary, the discussion for the short-time (unsteady) regime has been kept with some modifications, and the discussion for the long-time (steady) regime has been added as aforementioned.
Minor: the prefactor in (S15) is not correct, it should be 6π for a solid particle and less (depending on surface condition) for a bubble.
Thanks for the nice point! In our experiments, the Reynolds number for bubble motion = ≈ 10 2 , where the density of water = 10 3 kg/m 3 , the characteristic velocity for bubble motion =0.1 m/s, the characteristic length for bubble D=1 mm, and the viscosity of water = 1mPa • s. According to previous work on the viscous drag force on a spherical bubble [Physics of Fluids 10, 550 (1998)], when Re is large, the viscous contribution ( ) = 12 ( ) .
In the revised manuscript, the reference is added on page 8 as "the Stokes viscous drag force is raised with ( = −12 (31))". In the revised supplementary materials, the reference is added on page 9 as "the viscous force according to Ref. (5) is:"